Nonparametric Statistics
AY 2023-2024
Let us try to code the iterative algorithm described in slide 12, but with an additional (robust) scale estimate \(\hat{\sigma}\)
Consider the natural logs of the annual incomes of \(10\) people
manual_M_location <-
function(x,
k,
type = c("Huber", "Tukey"),
tol = sqrt(.Machine$double.eps),
itermax = 1000) {
mu <- mu_old <- mu_vec <- median(x) # initial value
rob_scale <- mad(x) # it will be kept fixed
crit <- TRUE
iter <- 0
weigth_f <- switch (type,
"Huber" = function(x) pmin(1, k/abs(x)),
"Tukey" = function(x) ((1-(x/k)^2)^2)*(abs(x)<=k)
)
while(crit){
w_i <- weigth_f((x-mu)/rob_scale)
mu <- weighted.mean(x = x,w = w_i)
err <- abs(mu-mu_old)
mu_old <- mu
mu_vec <- c(mu_vec,mu)
iter <- iter+1
crit <- (err > tol & iter < itermax)
}
list(mu=mu, s=rob_scale, w_i=w_i, it=iter)
}